###census api key
census_api_key("7492e696e8c0d9bcb8a075682c4dd9e68aad700a", install=TRUE, overwrite=TRUE)

###state  #population,income,white,hisp/lat
state <- get_acs(geography = "state", 
                       variables = c("B01003_001", #total pop
                                     "B19013_001", #med household income
                                     "B19301_001", #per capita income
                                     "B02001_002", #white alone (race)
                                     "B02001_001", #total pop (race)
                                     "B03003_003", #hisp/lat (hisp/lat)
                                     "B03003_001"), #total pop (hisp/lat)
                       year = 2019, 
                       geometry=TRUE,
                 output="wide") %>%
  filter(NAME!="Alaska" & NAME!="Hawaii" & NAME!="Puerto Rico") %>%
  rename(state_pop=B01003_001E,
         state_medHouseholdIncome=B19013_001E,
         state_perCapInc=B19301_001E,
         state_race_whitePop=B02001_002E,
         state_race_totalPop=B02001_001E,
         state_hisp_lat_yes=B03003_003E,
         state_hisp_lat_totalPop=B03003_001E) %>%
  mutate(state_prop_race_nonwhitePop=round((1-state_race_whitePop/state_race_totalPop),2),
         state_prop_hisp_lat_yes=round((state_hisp_lat_yes/state_hisp_lat_totalPop),2)) %>%
  st_transform(crs=4326) 

###relevant state variables to join to mines
state2 <- state %>%
  select(NAME,
         state_medHouseholdIncome,
         state_perCapInc,
         state_prop_race_nonwhitePop,
         state_prop_hisp_lat_yes) %>%
  st_drop_geometry()

###mines in U.S.  #cobalt, lithium, tellurium, rare earths
minesKeyMinerals <- minesRaw %>% #filter to key  minerals
  filter(grepl("Cobalt|Lithium|Tellurium|Rare-Earth Elements", CRITICAL_MINERAL)) %>%
  st_transform(crs=4326) %>%
  left_join(state2,by=c("LOC_DETAIL"="NAME"))
  
###create 35-mile buffer polygons around mines   #cobalt, lithium, tellurium, rare earths
buffer <- st_buffer(minesKeyMinerals, dist = 56327) %>% #in meters
  st_transform(crs=4326) 

###county  #population,income,white,hisp/lat
county <- get_acs(geography = "county",
                       variables = c("B01003_001", #total pop
                                     "B19013_001", #med household income
                                     "B19301_001", #per capita income
                                     "B02001_002", #white alone (race)
                                     "B02001_001", #total pop (race)
                                     "B03003_003", #hisp/lat
                                     "B03003_001"), #total pop (hisp/lat)
                       year = 2019,
                       geometry=TRUE,
                  output="wide") %>%
  separate(NAME,into=c("County","State"),sep=", ") %>%
  filter(State!="Alaska" & State!="Hawaii" & State!="Puerto Rico") %>%
  rename(pop=B01003_001E,
         medHouseholdIncome=B19013_001E,
         perCapInc=B19301_001E,
         race_whitePop=B02001_002E,
         race_totalPop=B02001_001E,
         hisp_lat_yes=B03003_003E,
         hisp_lat_totalPop=B03003_001E) %>% 
  mutate(prop_race_nonwhitePop=round((1-race_whitePop/race_totalPop),2),
         prop_hisp_lat_yes=round((hisp_lat_yes/hisp_lat_totalPop),2)) %>%
  st_transform(crs=4326)

###place   #population,income,white,hisp/lat
place <- get_acs(geography = "place", 
                       variables = c("B01003_001", #total pop
                                     "B19013_001", #med household income
                                     "B19301_001", #per capita income
                                     "B02001_002", #white alone (race)
                                     "B02001_001", #total pop (race)
                                     "B03003_003", #hisp/lat (hisp/lat)
                                     "B03003_001"), #total pop (hisp/lat)
                       year = 2019, 
                       geometry=TRUE,
                  output="wide") %>%
  separate(NAME,into=c("County","State"),sep=", ") %>%
  filter(State!="Alaska" & State!="Hawaii" & State!="Puerto Rico") %>%
  rename(pop=B01003_001E,
         medHouseholdIncome=B19013_001E,
         perCapInc=B19301_001E,
         race_whitePop=B02001_002E,
         race_totalPop=B02001_001E,
         hisp_lat_yes=B03003_003E,
         hisp_lat_totalPop=B03003_001E) %>% 
  mutate(prop_race_nonwhitePop=round((1-race_whitePop/race_totalPop),2),
         prop_hisp_lat_yes=round(hisp_lat_yes/hisp_lat_totalPop),2) %>%
  st_transform(crs=4326)

placeCentroid <- place %>% #places to centroids
  st_centroid()

intersectCentroid <- st_intersection(placeCentroid, buffer) #centroids in buffer

minesMapPlaceCentroids_buffer <- ggplot() +
  geom_polygon(data=base1, fill="white", color="grey45", size=0.3,
               aes(x=long, y=lat, group=group)) + #add base map
  geom_sf(data=intersectCentroid, color="grey70") + #add places
  geom_sf(data=buffer, fill=NA, linetype="dotted", color="red", size=0.3)+
  geom_point(data=minesKeyMinerals, color="slateblue", 
             aes(x=LONGITUDE, y=LATITUDE, 
                 text=paste(CRITICAL_MINERAL))) + #add mines
  theme_minimal()+ #remove grey background
  theme(legend.position="none", #remove legend, axes labels/ticks
        axis.title.x=element_blank(), 
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank()) +
  ggtitle("Select Places in 35-mile Buffer")

intersect <- st_intersection(county, buffer) %>%
  sf::st_cast("MULTIPOLYGON") #counties in buffer

fullCounties_buffer <- county %>%
  filter(GEOID %in% intersect$GEOID) #select full counties

mapMines_fullCounties_buffer <- ggplot() +
  geom_sf(data=county, color="grey45", size=0.3, fill=NA)+ #add base map
  geom_sf(data=fullCounties_buffer, color="grey45", size=0.3, fill="grey75")+ 
  geom_sf(data=buffer, fill=NA, linetype="dotted", size=0.3, color="red")+ #add buffer 
    geom_point(data=minesKeyMinerals, color="slateblue",
             aes(x=LONGITUDE, y=LATITUDE,
                 text=paste("",CRITICAL_MINERAL))) + #add mines
  theme_minimal()+ #remove grey background
  theme(legend.position="none", #remove legend, axes labels/ticks
        axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank()) +
  ggtitle("Select Counties in 35-mile Buffer")

mappingMines <- function(data,geom,bufferDist,variable,pointColor="tomato3") {
  
  ###create buffer polygons
  buffer <- st_buffer(data, dist = bufferDist) %>% #in meters
  st_transform(crs=4326) 
  
   ###filter to buffer
  if (geom=="place") {
    intersect <- st_intersection(placeCentroid, buffer)
  } else {
    intersect <- st_intersection(county, buffer)
  }
   
    groupedData <- intersect %>%
      group_by(CRITICAL_MINERAL,DEPOSIT_NAME,LONGITUDE,LATITUDE,LOC_DETAIL, 
               state_medHouseholdIncome,state_perCapInc,
               state_prop_race_nonwhitePop,state_prop_hisp_lat_yes) %>%
      summarize(population=sum(pop),
                medianHouseholdIncome=round(mean(medHouseholdIncome),2),
                perCapitaIncome=round(mean(perCapInc),2),
                propPopNonwhite=round(mean(prop_race_nonwhitePop),2),
                propHisLat=round(mean(prop_hisp_lat_yes),2)) %>%
      mutate(medHouseholdIncome_STATE=round(medianHouseholdIncome/state_medHouseholdIncome,2),
             perCapitaIncome_STATE= round(perCapitaIncome/state_perCapInc,2),
             propPopNonwhite_STATE=round(propPopNonwhite/state_prop_race_nonwhitePop,2),
            propHisLat_STATE=round(propHisLat/state_prop_hisp_lat_yes,2)) 
    
    ###filter to buffer
  if (geom=="place") {
    groupedData <- groupedData %>%
      sf::st_cast("POINT")
  } else {
    groupedData <- groupedData %>%
      sf::st_cast("MULTIPOLYGON")
  }
    
  ###set variables + labels for map
  if (variable=="population") {
    column <- groupedData$population
    result <- "population"
  } else if (variable=="median household income") {
    column <- groupedData$medianHouseholdIncome
    result <- "median household income"
  } else if (variable=="median household income2") {
    column <- groupedData$medHouseholdIncome_STATE
    result <- "med household inc (prop of state)"
  } else if (variable=="per capita income") {
    column <- groupedData$perCapitaIncome
    result <- "per capita income"
  } else if (variable=="per capita income2") {
    column <- groupedData$perCapitaIncome_STATE
    result <- "per capita income (prop of state)"
  } else if (variable=="race") {
    column <- groupedData$propPopNonwhite
    result <- "prop nonwhite pop"
  } else if (variable=="race2") {
    column <- groupedData$propPopNonwhite_STATE
    result <- "nonwhite pop (prop of state)"
  } else if (variable=="ethnicity") { 
    column <- groupedData$propHisLat
    result <- "prop Hisp/Lat pop"
  } else if (variable=="ethnicity2") { 
    column <- groupedData$propHisLat_STATE
    result <- "Hisp/Lat pop (prop of state)"
  }

  ###map
  map <- ggplot() +
    geom_polygon(data=base1, color = "grey18", fill="grey12", size=0.4,
                 aes(x=long, y=lat, group=group)) + #add base map
    geom_point(groupedData, color=pointColor, alpha=0.5,
            mapping=aes(x=LONGITUDE, y=LATITUDE, size=column,
                        label1=DEPOSIT_NAME, label2=CRITICAL_MINERAL, label3=LOC_DETAIL,
                        text=paste("name: ",DEPOSIT_NAME,"</br>",
                                   "</br>mineral: ",CRITICAL_MINERAL,
                                   "</br>state: ",LOC_DETAIL,
                                   "</br>",result,": ",column,sep=""))) +
    theme_minimal()+ #remove grey background
    theme(legend.position="none", #remove legend, axes labels/ticks
          axis.title.x=element_blank(),
          axis.text.x=element_blank(),
          axis.ticks.x=element_blank(),
          axis.title.y=element_blank(),
          axis.text.y=element_blank(),
          axis.ticks.y=element_blank()) +
    ggtitle(paste("Mines by",result))
  
  ggplotly(map,tooltip="text") %>%
    style(hoverinfo="none", traces=1)
  
}

mappingMines(data=minesKeyMinerals,geom="place",bufferDist=56327,variable="population")

intersectCentroid %>%
  as_tibble() %>%
  group_by(DEPOSIT_NAME,CRITICAL_MINERAL,LOC_DETAIL) %>%
  summarize(population=sum(pop)) %>%
  arrange(population)

mappingMines(data=minesKeyMinerals,geom="county",bufferDist=56327,
             variable="median household income",pointColor="darkseagreen4")

intersect %>%
  as_tibble() %>%
  group_by(DEPOSIT_NAME,CRITICAL_MINERAL,LOC_DETAIL) %>%
  summarize(medianHouseholdIncome=round(mean(medHouseholdIncome),2)) %>% #round to 2 decimal places
  arrange(medianHouseholdIncome)

intersect %>%
  as_tibble() %>%
  group_by(DEPOSIT_NAME,CRITICAL_MINERAL,LOC_DETAIL,state_medHouseholdIncome) %>%
  summarize(medianHouseholdIncome=round(mean(medHouseholdIncome),2)) %>%
  #round to 2 decimal places
  mutate(medHouseholdIncome_STATE=round(medianHouseholdIncome/state_medHouseholdIncome,2)) %>% 
  select(DEPOSIT_NAME,CRITICAL_MINERAL,LOC_DETAIL,medHouseholdIncome_STATE) %>% #remove extra variables
  arrange(medHouseholdIncome_STATE)

mappingMines(data=minesKeyMinerals,geom="county",bufferDist=56327,
             variable="race",pointColor="yellow4")

intersect %>%
  as_tibble() %>%
  group_by(DEPOSIT_NAME,CRITICAL_MINERAL,LOC_DETAIL) %>%
  summarize(propPopNonwhite=round(mean(prop_race_nonwhitePop),2)) %>% #round to 2 decimal places
  arrange(propPopNonwhite)

intersect %>%
  as_tibble() %>%
  group_by(DEPOSIT_NAME,CRITICAL_MINERAL,LOC_DETAIL,state_prop_race_nonwhitePop) %>%
  summarize(propPopNonwhite=round(mean(prop_race_nonwhitePop),2)) %>%
  mutate(propPopNonwhite_STATE=round(propPopNonwhite/state_prop_race_nonwhitePop,2)) %>% 
  #round to 2 decimal places
  select(DEPOSIT_NAME,CRITICAL_MINERAL,LOC_DETAIL,propPopNonwhite_STATE) %>% #remote extra variables
  arrange(propPopNonwhite_STATE)

mappingMines(data=minesKeyMinerals,geom="county",bufferDist=56327,
             variable="ethnicity",pointColor="darkorange4")

intersect %>%
  as_tibble() %>%
  group_by(DEPOSIT_NAME,CRITICAL_MINERAL,LOC_DETAIL) %>%
  summarize(propHisLat=round(mean(prop_hisp_lat_yes),2)) %>% #round to 2 decimal places
  arrange(propHisLat)

intersect %>%
  as_tibble() %>%
  group_by(DEPOSIT_NAME,CRITICAL_MINERAL,LOC_DETAIL,state_prop_hisp_lat_yes) %>%
  summarize(propHisLat=round(mean(prop_hisp_lat_yes),2)) %>% #round to 2 decimal places
  mutate(propHisLat_STATE=round(propHisLat/state_prop_hisp_lat_yes,2))%>%
  select(DEPOSIT_NAME,CRITICAL_MINERAL,LOC_DETAIL,propHisLat_STATE) %>%
  arrange(propHisLat_STATE)

base_caseSites <- state %>%
  filter(NAME %in% c("Nevada",
                     "Oregon",
                     "Idaho"))

minesKeyMinerals_caseSites <- minesKeyMinerals %>%
  filter(DEPOSIT_NAME %in% c("Kings Valley","Blackbird"))

buffer_caseSites <- buffer %>%
  filter(DEPOSIT_NAME %in% c("Kings Valley","Blackbird"))

intersectCentroid_caseSites <- st_intersection(placeCentroid, buffer_caseSites)  %>%#centroids in buffer
  mutate(medHouseholdIncome_STATE=round(medHouseholdIncome/state_medHouseholdIncome,2),
         perCapInc_STATE=round(perCapInc/state_perCapInc,2),
         propPopNonwhite_STATE=round(prop_race_nonwhitePop/state_prop_race_nonwhitePop,2),
         propHisLat_STATE=round(prop_hisp_lat_yes/state_prop_hisp_lat_yes,2)) %>%
  filter(pop>0)

placeCentroid_caseSites <- placeCentroid %>%
  filter(State %in% c("Nevada",
                      "Oregon",
                      "Idaho"))

minesMapPlaceCentroids_buffer_caseSites <- ggplot() +
  geom_sf(data=base_caseSites, fill="grey90", color="grey85", size=0.3) + #add base map
  geom_sf(data=placeCentroid_caseSites, color="grey85", fill=NA)+
  geom_sf(data=buffer_caseSites, fill=NA, linetype="dotted", color="salmon4", size=0.3)+
  geom_sf(data=intersectCentroid_caseSites, color="salmon3", alpha=0.5, 
          aes(size=pop,
          label1=medHouseholdIncome,label2=medHouseholdIncome_STATE,
          label3=perCapInc,label4=perCapInc_STATE,
          label5=prop_race_nonwhitePop,label6=propPopNonwhite_STATE,
          label7=prop_hisp_lat_yes,label8=propHisLat_STATE,
          label9=County,
          text=paste("place: ",County,"</br>",
                     "</br>population: ",pop,
                     "</br>income: ",medHouseholdIncome,
                     "</br>income (prop state): ",medHouseholdIncome_STATE,
                     "</br>per capita income: ",perCapInc,
                     "</br>per capita income (prop state): ",perCapInc_STATE,
                     "</br>prop nonwhite population: ",prop_race_nonwhitePop,
                     "</br>nonwhite population (prop state): ",propPopNonwhite_STATE,
                     "</br>prop Hisp/Lat: ",prop_hisp_lat_yes,
                     "</br>Hisp/Lat population (prop state): ",propHisLat_STATE,
                     sep=""))) + #add places
  geom_point(data=minesKeyMinerals_caseSites, color="grey30", 
             aes(x=LONGITUDE, y=LATITUDE,
                 label1=DEPOSIT_NAME,label2=CRITICAL_MINERAL,label3=LOC_DETAIL,
                 text=paste("name: ",DEPOSIT_NAME,"</br>",
                         "</br>mineral: ",CRITICAL_MINERAL,
                         "</br>state: ",LOC_DETAIL,sep=""))) + #add mines
  theme_minimal()+ #remove grey background
  theme(legend.position="none", #remove legend, axes labels/ticks
        axis.title.x=element_blank(), 
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()) +
  ggtitle("Places near Kings Valley and Blackbird deposits")

ggplotly(minesMapPlaceCentroids_buffer_caseSites,tooltip="text") %>%
    style(hoverinfo="none", traces=c(1,2))